function [TQ_tab,cwt_struct,tqplot] = ...
T_Q_stack(traces,... traces:irisFetch trace object
    Tmax,... %max period
    Tmin,... %min period
    minSTALTA,... %min amp/LTA ratio
    minpromstd,... %min prominence of peaks (rel to std of log LTA)
    minampstd,... %min amp of peaks (rel to std of log LTA)
    minTseprat,... %minimum ratio of periods between peaks
    mintsep,... %minimum separation between peaks in time
    fitcycles,... %sets number of cycles to exp fit
    minadjrsquare,... %min adjusted r squared value for exp fit
    VoicesPerOctave,... %for cwt
    WaveletWeights,... %[morse, bump, spec, narrowmorse], weights on different transforms
    QWaveletWeights,...%for expfit
    WaveletParameters,... %morse wavelet params for cwt
    NarrowWaveletParameters,...
    specwindowtime,... %spectrogram window 
    decaythresh,... %min remaining amplitude fraction per cycle
    Aweight,... %weight of mode "goodness" on amplitude
    gweight,... %weight mode "goodness" on decay rate
    rsquaredweight,... %weight mode "goodness" on exp fit r^2
    promweight,... %weight mode "goodness" on prominence
    startbounds,... %peak start time bounds
    makeplot,... %plot
    weight_decay_scale,... %decaying exp fit weights, smaller for faster decay
    c_quantile,... %quantile at which to set constant offset for exp fits
    c_std_frac,... %fraction of range beneath min to set c
    ampcycles,... %number of cycles to sum over to check that is a local max
    minhalflife,...%minumum oscillation half life 
    cwt_compress,... % output compressed cwt res {timewindow,numfreqbins} (empty for no compression)
    fitstartoffsetcycles,...%offset to start expfit,recommend 1-2 cycles
    LTAcycleoffset,... %amount to shift LTA (accouns for transform delay)
    LTA_window,... %windowspan for long-term avg
    expmodelvars) %variables to fit in "(d-c)*exp(bx)+c"
        %'bcd' for full 3 param or 
        %'bc' for fixed amp variable offset
        %'b' for fixed amp and offset
        %'bd' for fixed offset
        % 'b_under' for fixed amp + offset and fit under all data points


% calculate T and Q from stacked wavelet transforms
% Josh Crozier 2019

%some hardcoded parameters still scattered in code, could fix at some point


%% define fit models
if strcmp(expmodelvars,'bcd')
    %exp fit with all param variable
    expfitmodel = fittype('(d - c)*exp(b*x) + c',...
        'independent',{'x'},...
        'dependent',{'y'},...
        'coefficients',{'b','c','d'}); 
    
elseif strcmp(expmodelvars,'bc')
    %exp fit with fixed amp and variable offset
    expfitmodel = fittype('(d - c)*exp(b*x) + c',...
        'independent',{'x'},...
        'dependent',{'y'},...
        'coefficients',{'b','c'},...
        'problem',{'d'});
    
elseif strcmp(expmodelvars,'b')
    %exp fit with fixed amp and variable offset
    expfitmodel = fittype('(d - c)*exp(b*x) + c',...
        'independent',{'x'},...
        'dependent',{'y'},...
        'coefficients',{'b'},...
        'problem',{'c','d'});
    
elseif strcmp(expmodelvars,'bd')
    %will just use built in expfit
    
elseif strcmp(expmodelvars,'b_under')
    %just solve directly
    
else
    error('invalid exp fit model')
end

%% setup
if mod(traces(1).sampleCount,2) == 1
    error('odd data length')
end
Ts = seconds(1./traces(1).sampleRate);
time = traces(1).startTime + ...
    (0:traces(1).sampleCount-1).'*Ts;

startbounds = [time(1) + startbounds(1), time(end) - startbounds(2) + Ts];

min_g = log(0.5)/seconds(minhalflife);

%% make stacked cwt/spectrograms
count = 0;
specwindow = round(specwindowtime/seconds(Ts)/2)*2 + 1;
specnoverlap = specwindow-1;
for tracein = 1:length(traces) % loop over traces
    if traces(tracein).sampleCount ~= length(time)
        error('different data length')
    end
    try
        %% cwt transforms and spectrograms
        if count == 0
            %morse cwt
            if WaveletWeights(1) > 0 || QWaveletWeights(1) > 0
                [wt,s,~] = cwt(traces(tracein).data,'morse',1/seconds(Ts),...
                    'ExtendSignal',false,...
                    'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                    'VoicesPerOctave',VoicesPerOctave,...
                    'WaveletParameters',WaveletParameters);
                wt = (abs(wt));
            end
            %narrow morse
            if WaveletWeights(4) > 0 || QWaveletWeights(4) > 0
                [wtn,s,~] = cwt(traces(tracein).data,'morse',1/seconds(Ts),...
                    'ExtendSignal',false,...
                    'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                    'VoicesPerOctave',VoicesPerOctave,...
                    'WaveletParameters',NarrowWaveletParameters);
                wtn = (abs(wtn));
            end
            %bump cwt
            if WaveletWeights(2) > 0 || QWaveletWeights(2) > 0
                [wtb,s,~] = cwt(traces(tracein).data,'bump',1/seconds(Ts),...
                    'ExtendSignal',false,...
                    'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                    'VoicesPerOctave',VoicesPerOctave);
                wtb = (abs(wtb));
            end
            %spectrogram
            if ~exist('s','var')
                s = flipud(linspace(1/seconds(Tmax),1/seconds(Tmin),ceil(VoicesPerOctave*log2(Tmax/Tmin))).');
            end
            if WaveletWeights(3) > 0 || QWaveletWeights(3) > 0
                [sspec,~,timespec] = spectrogram(traces(tracein).data,...
                    specwindow,specnoverlap,s,1/seconds(Ts),'power');
                sspec = (abs(sspec));
            else
                timespec = time(time >= time(1) + specwindowtime/2 & ...
                    time <= time(end) - specwindowtime/2).';
            end
            
            T = seconds(1./s);

        else
            %morse cwt
            if WaveletWeights(1) > 0 || QWaveletWeights(1) > 0
                wt = wt + abs(cwt(traces(tracein).data,'morse',1/seconds(Ts),...
                    'ExtendSignal',false,...
                    'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                    'VoicesPerOctave',VoicesPerOctave,...
                    'WaveletParameters',WaveletParameters));
            end
            %narrow morse cwt
            if WaveletWeights(4) > 0 || QWaveletWeights(4) > 0
                wtn = wtn + abs(cwt(traces(tracein).data,'morse',1/seconds(Ts),...
                    'ExtendSignal',false,...
                    'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                    'VoicesPerOctave',VoicesPerOctave,...
                    'WaveletParameters',NarrowWaveletParameters));
            end
            % bump cwt
            if WaveletWeights(2) > 0 || QWaveletWeights(2) > 0
                wtb = wtb + abs(cwt(traces(tracein).data,'bump',1/seconds(Ts),...
                    'ExtendSignal',false,...
                    'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                    'VoicesPerOctave',VoicesPerOctave));
            end
            %spectrogram
            if WaveletWeights(3) > 0 || QWaveletWeights(3) > 0
                sspec = sspec + (abs(spectrogram(traces(tracein).data,...
                    specwindow,specnoverlap,s,1/seconds(Ts),'power')));
            end
            
        end %end if count == 0
        count = count + 1;
        
        %{
        %debug plots
        figure(1),clf
        %DFT
        subplot(2,3,1)
        [pxx,wxx] = periodogram(traces(tracein).data,[],[],1/Ts);
        pxx = pxx(wxx > 1/seconds(Tmax) & wxx < 1/seconds(Tmin));
        wxx = wxx(wxx > 1/seconds(Tmax) & wxx < 1/seconds(Tmin));
        loglog(1./wxx,pxx)
        xlim([Tmin,Tmax])
        ylim([min(pxx),max(pxx)])
        %waveform
        subplot(2,3,3)
        plot(time,traces(tracein).data)
        xlim([min(time),max(time)])
        ylim([min(traces(tracein).data),max(traces(tracein).data)])
        %cwt
        subplot(2,3,2)
        p = pcolor(time,T,log10(abs(wt_temp)));
        %p.FaceColor = 'interp';
        p.EdgeColor = 'none';
        colormap jet
        grid on
        hold on
        ylabel('Period')
        xlabel('Time (UTC)')
        set(gca,'yscale','log')
        grid on
        set(gca,'layer','top');
        colorbar
        %}
        
    catch ME
        disp(getReport(ME))
        disp(traces(1).startTime)
        disp(size(traces(tracein).data))
        disp(size(traces(1).data))
    end
end % end loop over traces


%% combine different transforms
%remove cone of edge influence from cwts
% for jj = 1:size(wt,2)
%     wt(s<coi(jj),jj) = min(wt(:));
%     wtb(s<coib(jj),jj) = min(wtb(:));
% end

%trim to spectrogram edges
%!!!!debug catch
try
    timespecin = time >= timespec(1) & time <= timespec(end);
catch ME
    disp(size(time))
    disp(size(timespec))
    disp([time(1),time(end)])
    disp([timespec(1),timespec(end)])
    rethrow(ME)
end
if ~makeplot
    clear time
end

if sum(timespecin) > length(timespec)
    timespecin(find(timespecind, 1)) = false;
elseif sum(timespecin) < length(timespec)
    timespecin(find(timespecin, 1)-1) = true;
end

%combine
comb = zeros(length(s),length(timespec));
combQ = zeros(length(s),length(timespec));
if WaveletWeights(1) > 0 || QWaveletWeights(1) > 0
    wt = wt(:,timespecin);
    if WaveletWeights(1) > 0
        comb = comb + (wt/std(wt(:)))*WaveletWeights(1);
    end
    if QWaveletWeights(1) > 0
        combQ = combQ + (wt/std(wt(:)))*QWaveletWeights(1);
    end
    wt = log10(wt);
end
if WaveletWeights(4) > 0 || QWaveletWeights(4) > 0
    wtn = wtn(:,timespecin);
    if WaveletWeights(4) > 0
        comb = comb + (wtn/std(wtn(:)))*WaveletWeights(4);
    end
    if QWaveletWeights(4) > 0
        combQ = combQ + (wtn/std(wtn(:)))*QWaveletWeights(4);
    end
    wtn = log10(wtn);
end
if WaveletWeights(2) > 0 || QWaveletWeights(2) > 0
    wtb = wtb(:,timespecin);
    if WaveletWeights(2) > 0
        comb = comb + (wtb/std(wtb(:)))*WaveletWeights(2);
    end
    if QWaveletWeights(2) > 0
        combQ = combQ + (wtb/std(wtb(:)))*QWaveletWeights(2);
    end
    wtb = log10(wtb);
end
if WaveletWeights(3) > 0 || QWaveletWeights(3) > 0
    if WaveletWeights(3) > 0
        comb = comb + (sspec/std(sspec(:)))*WaveletWeights(3);
    end
    if QWaveletWeights(3) > 0
        combQ = combQ + (sspec/std(sspec(:)))*QWaveletWeights(3);
    end
    sspec = log10(sspec);
end
comb = log10(comb);
combQ = log10(combQ);
if ~makeplot
    clear wt
    clear wtb
    clear wtn
    clear sspec
end

%% find possible peaks
startminind = find(timespec>=startbounds(1),1);
if isempty(startminind)
    startminind = 1;
end
startmaxind = find(timespec>=startbounds(2),1);
if isempty(startmaxind)
    startmaxind = length(timespec);
end

% old method using whole window
%{
% stdamp = nanstd((comb(:,1:noiseind)),0,2);
% stdamp = movmax(stdamp,min_amp_std_T_span);
% comb_bounded = comb(:,startminind:startmaxind);
% minprom = stdamp*minpromstd;
% minprom = repmat(minprom,1,size(comb_bounded,2));

% %     Amin = max(wt_trim,2)*Aminscale;
% % [peakind1,prom1] = islocalmax(wt_trim(:,startminind:startmaxind),1,...
% %     'MinSeparation',0, 'MinProminence',0);
% [peakind2,prom2] = islocalmax(comb_bounded,2,...
%     'MinSeparation',floor(mintsep/Ts), 'MinProminence',0);
% peakind = peakind2 & prom2 >= minprom;
%
% %make table
% [Tind_temp,startind_temp] = find(peakind);
% peaklinind = sub2ind(size(comb_bounded),Tind_temp,startind_temp);
% TQ_tab = table(T(Tind_temp),'VariableNames',{'T'});
% TQ_tab.stdamp = stdamp(Tind_temp);
% TQ_tab.prom = prom2(peaklinind);
% TQ_tab.Tind = Tind_temp;
% TQ_tab.startind = startind_temp + startminind - 1;
% TQ_tab.start = timespec(TQ_tab.startind).';
% TQ_tab.Acwt = comb_bounded(peaklinind);
% 
% %remove small amp
% medamp = nanmedian(comb(:,1:noiseind),2);
% medamp = movmax(medamp,min_amp_std_T_span);
% TQ_tab.medamp = medamp(TQ_tab.Tind);
% minamp_samples = TQ_tab.medamp + TQ_tab.stdamp*minampstd;
% keepind = TQ_tab.Acwt >= minamp_samples;
% TQ_tab = TQ_tab(keepind,:);
% 
% %store for later use
% meanamp = nanmean(comb(:,1:noiseind),2);
% meanamp = movmax(meanamp,min_amp_std_T_span);
% TQ_tab.meanamp = meanamp(TQ_tab.Tind);
% minamp = nanmin(comb,[],2);
% minamp = movmin(minamp,min_amp_std_T_span);
% TQ_tab.minamp = minamp(TQ_tab.Tind);
% maxamp = nanmax(comb,2);
% maxamp = movmax(maxamp,min_amp_std_T_span);
% TQ_tab.maxamp = maxamp(TQ_tab.Tind);
% TQ_tab.rangeamp = TQ_tab.maxamp-TQ_tab.minamp;
%}


%moving mean and std
LTA = movmean(comb,[ceil(seconds(LTA_window/seconds(Ts))),0],2);
LT_std = movstd(comb,[ceil(seconds(LTA_window/seconds(Ts))),0],0,2);
LT_med = movmedian(comb,[ceil(seconds(LTA_window/seconds(Ts))),0],2);
LT_min = movmin(comb,[ceil(seconds(LTA_window/seconds(Ts))),0],2);
LT_max = movmax(comb,[ceil(seconds(LTA_window/seconds(Ts))),0],2);
%shift moving mean due to transform delay
for Tin = 1:length(T)
    offset = floor(seconds(T(Tin))*LTAcycleoffset/seconds(Ts));
    LTA(Tin,:) = [zeros(1,offset)+LTA(Tin,1), LTA(Tin,1:end-offset)];
    LT_std(Tin,:) = [zeros(1,offset)+LT_std(Tin,1), LT_std(Tin,1:end-offset)];
    LT_med(Tin,:) = [zeros(1,offset)+LT_med(Tin,1), LT_med(Tin,1:end-offset)];
    LT_min(Tin,:) = [zeros(1,offset)+LT_min(Tin,1), LT_min(Tin,1:end-offset)];
    LT_max(Tin,:) = [zeros(1,offset)+LT_max(Tin,1), LT_max(Tin,1:end-offset)];
end
%find possible peaks
[peakind2,prom2] = islocalmax(comb,2,...
    'MinSeparation',floor(seconds(mintsep/seconds(Ts))), 'MinProminence',0);
A_above_LTA = comb - LTA;
A_above_LTA(:,1:ceil(seconds(LTA_window/seconds(Ts)))) = -inf;
A_above_LTA(:,1:startminind) = -inf;
A_above_LTA(:,startmaxind:end) = -inf;
peakind = peakind2 & ...
    prom2 >= LT_std*minpromstd & ...
    A_above_LTA > LT_std*minampstd &...
    comb - LTA > log10(minSTALTA);
clear peakind2
if ~makeplot
    clear A_above_LTA
end

%make table
[Tind_temp,startind_temp] = find(peakind);
clear peakind
peaklinind = sub2ind(size(LTA),Tind_temp,startind_temp);
TQ_tab = table(T(Tind_temp),'VariableNames',{'T'});
LT_std = LT_std(peaklinind);
TQ_tab.stdamp = LT_std;
LTA = LTA(peaklinind);
TQ_tab.meanamp = LTA;
LT_med = LT_med(peaklinind);
TQ_tab.medamp = LT_med;
LT_min = LT_min(peaklinind);
TQ_tab.minamp = LT_min;
LT_max = LT_max(peaklinind);
TQ_tab.maxamp = LT_max;
TQ_tab.rangeamp = TQ_tab.maxamp-TQ_tab.minamp;
TQ_tab.prom = prom2(peaklinind);
clear prom2
TQ_tab.Tind = Tind_temp;
TQ_tab.startind = startind_temp;
TQ_tab.start = timespec(TQ_tab.startind).';
TQ_tab.Acwt = comb(peaklinind);

% if any([TQ_tab.start] < startbounds(1)) || any([TQ_tab.start] > startbounds(2))
%     'debug pause'
% end


%remove temporally close peaks
if makeplot
    TQ_tab_00 = TQ_tab;
end
keepind = true(size(TQ_tab.T));
for kk = 1:length(TQ_tab.T)
    for jj = kk+1:length(TQ_tab.T)
        if TQ_tab.Tind(kk) == TQ_tab.Tind(jj) ...
                && abs(TQ_tab.start(kk) - TQ_tab.start(jj)) < mintsep
            if TQ_tab.Acwt(kk) >= TQ_tab.Acwt(jj) && keepind(kk)
                keepind(jj) = false;
            else
                keepind(kk) = false;
            end
        end
    end
end
TQ_tab = TQ_tab(keepind,:);

%remove close peaks that aren't local max in freq over some timespan
if makeplot
    TQ_tab_0 = TQ_tab;
end
keepind = true(size(TQ_tab.T));
for kk = 1:length(TQ_tab.T)
    sumspan = TQ_tab.T(kk)*ampcycles;
    sumspan = ceil(seconds(sumspan/seconds(Ts)));
    %periods within some range of each potential peak
    T_upper = minTseprat*TQ_tab.T(kk);
    T_upper_in = discretize(T_upper, ...
        [seconds(0);(T(1:end-1) + T(2:end))/2;Inf]);
    T_lower = TQ_tab.T(kk)/minTseprat;
    T_lower_in = discretize(T_lower, [0;(T(1:end-1) + T(2:end))/2;seconds(Inf)]);
    total_amp_vec = sum(10.^(comb(T_lower_in:T_upper_in,...
        TQ_tab.startind(kk):TQ_tab.startind(kk) + sumspan)), 2);
    rel_in = TQ_tab.Tind(kk) - T_lower_in + 1;
    
    if any(total_amp_vec > total_amp_vec(rel_in))
        keepind(kk) = false;
    end
end
TQ_tab = TQ_tab(keepind,:);

%% calculate Q by exponetial regression
TQ_tab.Q = NaN(size(TQ_tab.T));
TQ_tab.g = NaN(size(TQ_tab.T));
TQ_tab.Afit = NaN(size(TQ_tab.T));
TQ_tab.cfit = NaN(size(TQ_tab.T));
%store goodness parameters
TQ_tab.rsquare = NaN(size(TQ_tab.T));
TQ_tab.dfe = NaN(size(TQ_tab.T));
TQ_tab.adjrsquare = NaN(size(TQ_tab.T));
TQ_tab.rmse = NaN(size(TQ_tab.T));

keepind = true(size(TQ_tab.T));
for kk = 1:length(TQ_tab.T)
    %start after peak due to smoothing effects of transforms
    startoffset = floor(fitstartoffsetcycles*seconds(TQ_tab.T(kk))/seconds(Ts));
    endin = min([floor(seconds(TQ_tab.T(kk))*fitcycles/seconds(Ts)) + TQ_tab.startind(kk),...
        size(combQ, 2)]);
    wtlocal = 10.^(combQ(TQ_tab.Tind(kk),TQ_tab.startind(kk) + startoffset:endin));
    %wtlocal = wtlocal(~isnan(wtlocal)).';
    tlocal = (seconds(0):Ts:Ts*(length(wtlocal)-1));  
    Weights = exp(-seconds(tlocal)/(weight_decay_scale*fitcycles*seconds(TQ_tab.T(kk))));    
    
    if any(wtlocal > wtlocal(1))
        keepind(kk) = false;
        TQ_tab.Q(kk) = NaN;
        TQ_tab.g(kk) = NaN;
        TQ_tab.cfit(kk) = NaN;
        TQ_tab.Afit(kk) = NaN;
        TQ_tab.rsquare(kk) = NaN;
        TQ_tab.dfe(kk) = NaN;
        TQ_tab.adjrsquare(kk) = NaN;
        TQ_tab.rmse(kk) = NaN;
    else
        if strcmp(expmodelvars,'bcd')
            %exp fit with offset
            scale = range(10.^(combQ(TQ_tab.Tind(kk),:)));
            [fitlocal,goflocal,~] = fit(seconds(tlocal)',wtlocal',expfitmodel,...
                'StartPoint',[0.25*log(decaythresh)/seconds(TQ_tab.T(kk)), ...
                    10.^(quantile((combQ(TQ_tab.Tind(kk),:)),c_quantile)),...
                    0.01*range(wtlocal) + wtlocal(1)],...
                'Lower',[2*log(decaythresh)/seconds(TQ_tab.T(kk)), ...
                    10.^(quantile((combQ(TQ_tab.Tind(kk), :)), 0.05)),...
                    wtlocal(1)],...
                'Upper',[0,...
                    10.^(quantile((combQ(TQ_tab.Tind(kk),:)),0.5)),...
                    0.1*range(wtlocal) + wtlocal(1)],...
                'Weights',Weights,...
                'MaxFunEvals',1000,...
                'MaxIter',800,...
                'Robust','off',...
                'DiffMaxChange',0.1,...
                'DiffMinChange',1E-9*scale,...
                'TolFun',1E-11*scale,...
                'TolX',1E-9*scale,...
                'Normalize','off');
            coeffvals = coeffvalues(fitlocal);
            TQ_tab.g(kk) = coeffvals(1);
            TQ_tab.cfit(kk) = coeffvals(2);
            TQ_tab.Afit(kk) = coeffvals(3);
            TQ_tab.rsquare(kk) = goflocal.rsquare;
            TQ_tab.dfe(kk) = goflocal.dfe;
            TQ_tab.adjrsquare(kk) = goflocal.adjrsquare;
            TQ_tab.rmse(kk) = goflocal.rmse;

        elseif strcmp(expmodelvars,'bc')
            %exp fit with fixed amp and variable offset
            scale = range(10.^(combQ(TQ_tab.Tind(kk),:)));
            [fitlocal,goflocal,~] = fit(seconds(tlocal)',wtlocal',expfitmodel,...
                'StartPoint',[0.5*log(decaythresh)/seconds(TQ_tab.T(kk)), ...
                    10.^(quantile((combQ(TQ_tab.Tind(kk),:)),c_quantile))],...
                'Lower',[2*log(decaythresh)/seconds(TQ_tab.T(kk)), ...
                    10.^(quantile((combQ(TQ_tab.Tind(kk), :)), 0.05))],...
                'Upper',[0,...
                    10.^(quantile((combQ(TQ_tab.Tind(kk),:)),0.5))],...
                'problem',wtlocal(1),...
                'Weights',Weights,...
                'MaxFunEvals',1000,...
                'MaxIter',800,...
                'Robust','off',...
                'DiffMaxChange',0.1,...
                'DiffMinChange',1E-9*scale,...
                'TolFun',1E-11*scale,...
                'TolX',1E-9*scale,...
                'Normalize','off');
            coeffvals = coeffvalues(fitlocal);
            TQ_tab.g(kk) = coeffvals(1);
            TQ_tab.cfit(kk) = coeffvals(2);
            TQ_tab.rsquare(kk) = goflocal.rsquare;
            TQ_tab.dfe(kk) = goflocal.dfe;
            TQ_tab.adjrsquare(kk) = goflocal.adjrsquare;
            TQ_tab.rmse(kk) = goflocal.rmse;

        elseif strcmp(expmodelvars,'b')
            %exp fit with fixed amp and  offset
            if isinf(c_quantile)
                c = 0;
            else
                c = 10.^(quantile((combQ(TQ_tab.Tind(kk),:)),c_quantile)); %could adjust
            end
            scale = range(10.^(combQ(TQ_tab.Tind(kk),:)));
            [fitlocal,goflocal,~] = fit(seconds(tlocal)',wtlocal',expfitmodel,...
                'StartPoint',0.5*log(decaythresh)/seconds(TQ_tab.T(kk)),...
                'Lower',2*log(decaythresh)/seconds(TQ_tab.T(kk)),...
                'Upper',0,...
                'problem',{c,...
                    wtlocal(1)},...
                'Weights',Weights,...
                'MaxFunEvals',1000,...
                'MaxIter',800,...
                'Robust','off',...
                'DiffMaxChange',0.1,...
                'DiffMinChange',1E-9*scale,...
                'TolFun',1E-11*scale,...
                'TolX',1E-9*scale,...
                'Normalize','off');
            coeffvals = coeffvalues(fitlocal);
            TQ_tab.g(kk) = coeffvals(1);
            TQ_tab.rsquare(kk) = goflocal.rsquare;
            TQ_tab.dfe(kk) = goflocal.dfe;
            TQ_tab.adjrsquare(kk) = goflocal.adjrsquare;
            TQ_tab.rmse(kk) = goflocal.rmse;

        elseif strcmp(expmodelvars,'bd')
            %exp fit with fixed offset
            if isinf(c_quantile)
                c = 0;
            else
                c = 10.^(quantile((combQ(TQ_tab.Tind(kk),:)),c_quantile)); %could adjust
            end
            scale = range(10.^(combQ(TQ_tab.Tind(kk),:)));
    %         [fitlocal,goflocal] = fit(tlocal,wtlocal - c,'exp1',...
    %             'StartPoint',[0.01*range(wtlocal) + wtlocal(1) - c,...
    %                 0.25*log(decaythresh)/TQ_tab.T(kk)],...
    %             'Lower',[wtlocal(1) - c,...
    %                 2*log(decaythresh)/TQ_tab.T(kk)],...
    %             'Upper',[0.1*range(wtlocal) + wtlocal(1) - c,...
    %                 0],...
    %             'Weights',Weights,...
    %             'MaxFunEvals',1000,...
    %             'MaxIter',800);
            [fitlocal,goflocal,~] = fit(seconds(tlocal)', wtlocal' - c,'exp1',...
                'StartPoint',[0.01*range(wtlocal) + wtlocal(1) - c,...
                    1*log(decaythresh)/seconds(TQ_tab.T(kk))],...
                'Lower',[wtlocal(1) - c,...
                    2*log(decaythresh)/seconds(TQ_tab.T(kk))],...
                'Upper',[0.1*range(wtlocal) + wtlocal(1) - c,...
                    -1E-6],...
                'Weights',Weights,...
                'MaxFunEvals',1000,...
                'MaxIter',800,...
                'Robust','off',...
                'DiffMaxChange',0.1,...
                'DiffMinChange',1E-9*scale,...
                'TolFun',1E-11*scale,...
                'TolX',1E-9*scale,...
                'Normalize','off');
            coeffvals = coeffvalues(fitlocal);
            TQ_tab.g(kk) = coeffvals(1);
            TQ_tab.Afit(kk) = coeffvals(2);
            TQ_tab.rsquare(kk) = goflocal.rsquare;
            TQ_tab.dfe(kk) = goflocal.dfe;
            TQ_tab.adjrsquare(kk) = goflocal.adjrsquare;
            TQ_tab.rmse(kk) = goflocal.rmse;

        elseif strcmp(expmodelvars,'b_under')
            %solve for highest b (slowest decay) that is under all data
            if isinf(c_std_frac)
                c = 0;
            else
                c = 10.^(TQ_tab.minamp(kk)-c_std_frac*TQ_tab.rangeamp(kk));
            end
            d = wtlocal(1);
            dataind = wtlocal<wtlocal(1);
            %datalen = sum(dataind);
            g = log((wtlocal(dataind) - c)/(d - c))./seconds(tlocal(dataind));
            g = nanmin(g);
            if wtlocal(2) > wtlocal(1)
                'bad peak'
            end
            g = g(1); %in case non-unique min
            TQ_tab.g(kk) = g;
            if TQ_tab.g(kk) > 0
                TQ_tab.g(kk) = 0;
            end

            wtpred = (d - c)*exp(g*seconds(tlocal)) + c;
            sse = sum((wtlocal - wtpred).^2);
            sst = (length(wtlocal)-2)*var(wtlocal);
            p = 1;%number model params
            TQ_tab.rsquare(kk) = 1 - sse/sst;
            TQ_tab.adjrsquare(kk) = 1 - sse/sst*(length(wtlocal)-2)/(length(wtlocal)-1-p);
            TQ_tab.rmse(kk) = sqrt((sse)/(length(wtlocal)-1));
            TQ_tab.dfe(kk) = length(wtlocal)-1-p;



            %'debug plot'
            %{
            if TQ_tab.T(kk) > seconds(2) %seconds(30) && TQ_tab.T(kk) < seconds(50)
                %'debug plot'
                figure(7),clf
                plot_startin = max([TQ_tab.startind(kk)-2*startoffset, 1]);
                plot_endin = min([floor(TQ_tab.T(kk)*fitcycles/Ts)+...
                    TQ_tab.startind(kk)+8*startoffset, size(combQ, 2)]);
                wtlocal_plot = 10.^combQ(TQ_tab.Tind(kk), plot_startin:plot_endin);
                tlocal_plot = [-fliplr(tlocal(2:startoffset+1+TQ_tab.startind(kk)-plot_startin)), ...
                    tlocal, (1:plot_endin-endin)*Ts + tlocal(end)];
                wtpred_plot = (d - c)*exp(g*seconds(tlocal_plot)) + c;
                plot(tlocal_plot,wtlocal_plot,'-b')
                hold on
                plot(tlocal,wtlocal,'b-','LineWidth',2)
                plot(tlocal_plot,wtpred_plot,'r-')
                Q = -pi/(seconds(TQ_tab.T(kk))*TQ_tab.g(kk));
                title(strcat('T=',num2str(seconds(TQ_tab.T(kk))),' Q=',num2str(Q)))
                grid on
                xlabel('time (s)')
                ylabel('Scalogram Amplitude')
%                 Qextra = 25;
%                 gextra = -pi/Qextra/seconds(TQ_tab.T(kk));
%                 wtpred_plot_extra = (d - c)*exp(gextra*seconds(tlocal_plot)) + c;
%                 Qextra2 = 35;
%                 gextra2 = -pi/Qextra2/seconds(TQ_tab.T(kk));
%                 wtpred_plot_extra2 = (d - c)*exp(gextra2*seconds(tlocal_plot)) + c;
%                 hold on
%                 plot(tlocal_plot,wtpred_plot_extra,'g')
%                 plot(tlocal_plot,wtpred_plot_extra2,'c')
                legend('Scalogram Amplitude','Section being fit',strcat('Q =',num2str(Q)));%,...
                    %strcat('Q =',num2str(Qextra)),strcat('Q =',num2str(Qextra2)))
                hold on
            end
            %}

        else
            error('invaled exp fit model')
        end %end choose fit model

        %!!!!debug plot!!!!
        %{
        if TQ_tab.T(kk) > seconds(30) && TQ_tab.T(kk) < seconds(50)
            figure(8),clf
            subplot(2,1,1)
            coeffvals = coeffvalues(fitlocal);
            plot_endin = min([floor(TQ_tab.T(kk)*fitcycles/Ts)+...
                TQ_tab.startind(kk)+4*startoffset, size(combQ, 2)]);
            wtlocal_plot = 10.^combQ(TQ_tab.Tind(kk), TQ_tab.startind(kk):plot_endin);
            tlocal_plot = [-fliplr(tlocal(2:startoffset+1)), ...
                tlocal, (1:plot_endin-endin)*Ts + tlocal(end)];
            if strcmp(expmodelvars,'b')
                plot(fitlocal,seconds(tlocal_plot),wtlocal_plot),hold on
                plot(seconds(tlocal),wtlocal,'b','LineWidth',2)
                title(strcat('T= ',num2str(seconds(TQ_tab.T(kk))),...
                    ' Q=',num2str(-pi/seconds(TQ_tab.T(kk)*coeffvals(1)))))
            else
                plot(fitlocal,seconds(tlocal_plot),wtlocal_plot),hold on
                plot(seconds(tlocal),wtlocal,'b','LineWidth',2)
                title(strcat('T= ',num2str(seconds(TQ_tab.T(kk))),...
                    ' Q=',num2str(-pi/seconds(TQ_tab.T(kk)*coeffvals(2)))))
            end
            subplot(2,1,2)
            plot(tlocal,Weights)
            title('weights')
            xlabel('t')
            ylabel('weights')
        end
        %}

        TQ_tab.Q(kk) = -pi/(seconds(TQ_tab.T(kk))*TQ_tab.g(kk));

        if TQ_tab.adjrsquare < minadjrsquare
            keepind(kk) = false;
        elseif TQ_tab.g(kk) < log(decaythresh)/seconds(TQ_tab.T(kk))
            keepind(kk) = false;
        elseif TQ_tab.g(kk) < min_g %not enough oscillations
            keepind(kk) = false;
        end
        
    end %end if all points not below first
end %end exp reg
temp_r2 = TQ_tab.adjrsquare + 1;
temp_r2(temp_r2 < 0) = 0;
% TQ_tab.goodness = (TQ_tab.Acwt+nanmin(wt_trim(:))).^(Aweight) ...
%     .*(exp(TQ_tab.g*TQ_tab.T(kk))).^(gweight) ...
%     .*(temp_r2).^(rsquaredweight) ...
%     .*TQ_tab.prom.^(promweight);
TQ_tab.goodness = (exp(TQ_tab.Acwt)).^(Aweight) ...
    .*(exp(TQ_tab.g.*seconds(TQ_tab.T))).^(gweight) ...
    .*(temp_r2).^(rsquaredweight) ...
    .*TQ_tab.prom.^(promweight);


if makeplot
    TQ_tab_1 = TQ_tab;
end

%remove bad peaks
%Tind = Tind(keepind);
%startind = startind(keepind);
TQ_tab = TQ_tab(keepind,:);


%% sort through possible peaks
%remove close peaks in time or freq
if makeplot
    TQ_tab_2 = TQ_tab;
end
keepind = true(size(TQ_tab.T));
for kk = 1:length(TQ_tab.T)
    for jj = kk+1:length(TQ_tab.T)
        if abs(TQ_tab.start(kk) - TQ_tab.start(jj)) < mintsep
            if (TQ_tab.T(kk) >= TQ_tab.T(jj) && TQ_tab.T(kk)/TQ_tab.T(jj) < minTseprat) ||...
                    (TQ_tab.T(kk) < TQ_tab.T(jj) && TQ_tab.T(jj)/TQ_tab.T(kk) < minTseprat) 
                if TQ_tab.goodness(kk) > TQ_tab.goodness(jj)
                    keepind(jj) = false;
                elseif TQ_tab.goodness(kk) < TQ_tab.goodness(jj)
                    keepind(kk) = false;
                else
                    keepind(jj) = false; %arbitrary
                end
            end
        end
    end
end
%}
%remove boundary peaks (often artifacts)
keepind(TQ_tab.T == T(end) | TQ_tab.T == T(1)) = false;
%     Tind = Tind(keepind);
%     startind = startind(keepind);
TQ_tab = TQ_tab(keepind,:);

%store startTimes
TQ_tab.windowstartTime = repmat(traces(1).startTime,height(TQ_tab),1);
TQ_tab.channelCount = repmat(length(traces),height(TQ_tab),1);


%% plots
if makeplot
    if sum(WaveletWeights>0) == 1
        plotrows = 2;
        plotcols = 1;
    elseif sum(WaveletWeights>0) < 4
        plotrows = 2;
        plotcols = 1;
    else
        plotrows = 3;
        plotcols = 2;
    end
    plotcount = 0;
    tqplot = figure('name','TQwavelet');
%     tplot = timespec;%time range to plot
%     minutein = find(tplot.Minute ~= tplot(1).Minute & ...
%         mod(tplot.Minute,2) == 0, 1);
%     xticklist = minutein-120:60:60*60*4;
%     xticklabellist = datestr(traces(1).startTime + seconds(xticklist),'HH:MM');
%     xticklabellist(2:5:end,:) = ' ';
%     xticklabellist(3:5:end,:) = ' ';
%     xticklabellist(4:5:end,:) = ' ';
%     xticklabellist(5:5:end,:) = ' ';
    %% cwt morse
    %
    if WaveletWeights(1) > 0
        plotcount = plotcount + 1;
        subplot(plotrows,plotcols,plotcount)
        p = pcolor(timespec,seconds(T),(wt));
        caxis([quantile((wt(:)),0.05),quantile((wt(:)),1)])
        %p.FaceColor = 'interp';
        p.EdgeColor = 'none';
        colormap jet
        grid on
        hold on
    %     coip = fill([time(1);time;time(end);time(1)],1./[min(s);coi;min(s);min(s)],'w');
    %     coip.FaceAlpha = 0.8;
        ylabel('Period (s)')
        xlabel('Time (UTC)')
        set(gca,'yscale','log')
        grid on
        set(gca,'layer','top');
        %xticks(xticklist)
        xtickangle(0)%, grid minor
        %xticklabels(xticklabellist)
        title(strcat('Scalogram'))
        yticks([1:1:5,10:10:50])
        
        if ~isempty(TQ_tab)
%             scatter(TQ_tab_00.start,seconds(TQ_tab_00.T),16,'w','.')
%             scatter(TQ_tab_0.start,seconds(TQ_tab_0.T),16,'w','*')
%             scatter(TQ_tab_1.start,seconds(TQ_tab_1.T),36,'m','.')
%             scatter(TQ_tab_2.start,seconds(TQ_tab_2.T),36,'m','*')
            scatter(TQ_tab.start,seconds(TQ_tab.T),36,'k','o')
    %         text(TQ_tab_1.start+10,TQ_tab_1.T,strcat(...
    %             'Q=',num2str(TQ_tab_1.Q,'%.1f'),...
    %         '  p=',num2str(TQ_tab_1.goodness,'%.2e')),'Color','w')
    %         text(TQ_tab_2.start+10,TQ_tab_2.T,strcat(...
    %             'Q=',num2str(TQ_tab_2.Q,'%.1f'),...
    %         '  p=',num2str(TQ_tab_2.goodness,'%.2e')),'Color','m')
            text(TQ_tab.start+seconds(30),seconds(TQ_tab.T),strcat(...
                'Q=',num2str(TQ_tab.Q,'%.2f')),'Color','k','FontWeight','bold')
        end
    end
    %}
    
    %% cwt narrow morse
    %{
    if WaveletWeights(4) > 0
        plotcount = plotcount + 1;
        subplot(plotrows,plotcols,plotcount)
        p = pcolor(timespec,seconds(T),(wtn));
        caxis([quantile((wtn(:)),0.05),quantile((wtn(:)),1)])
        %p.FaceColor = 'interp';
        p.EdgeColor = 'none';
        colormap jet
        grid on
        hold on
    %     coip = fill([time(1);time;time(end);time(1)],1./[min(s);coi;min(s);min(s)],'w');
    %     coip.FaceAlpha = 0.8;
        ylabel('Period (s)')
        xlabel('Time (UTC)')
        set(gca,'yscale','log')
        grid on
        set(gca,'layer','top');
        %xticks(xticklist)
        xtickangle(0)%, grid minor
        %xticklabels(xticklabellist)
        title(strcat('Scalogram (narrow)'))
        yticks([1:1:5,10:10:50])
        
        if ~isempty(TQ_tab)
%             scatter(TQ_tab_00.start,seconds(TQ_tab_00.T),16,'w','.')
%             scatter(TQ_tab_0.start,seconds(TQ_tab_0.T),16,'w','*')
%             scatter(TQ_tab_1.start,seconds(TQ_tab_1.T),36,'m','.')
%             scatter(TQ_tab_2.start,seconds(TQ_tab_2.T),36,'m','*')
            scatter(TQ_tab.start,seconds(TQ_tab.T),36,'k','o')
    %         text(TQ_tab_1.start+10,TQ_tab_1.T,strcat(...
    %             'Q=',num2str(TQ_tab_1.Q,'%.1f'),...
    %         '  p=',num2str(TQ_tab_1.goodness,'%.2e')),'Color','w')
    %         text(TQ_tab_2.start+10,TQ_tab_2.T,strcat(...
    %             'Q=',num2str(TQ_tab_2.Q,'%.1f'),...
    %         '  p=',num2str(TQ_tab_2.goodness,'%.2e')),'Color','m')
            text(TQ_tab.start+seconds(30),seconds(TQ_tab.T),strcat(...
                'Q=',num2str(TQ_tab.Q,'%.2f')),'Color','k','FontWeight','bold')
        end
    end
    %}
    
    %% cwt bump
    if WaveletWeights(2) > 0
        plotcount = plotcount + 1;
        subplot(plotrows,plotcols,plotcount)
%     count = 0;
%     for tracein = 1:length(traces) % loop over traces
%         if count == 0
%             [wtb_temp,s,coib] = cwt(traces(tracein).data,'bump',1/Ts,...
%                 'ExtendSignal',false,...
%                 'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
%                 'VoicesPerOctave',VoicesPerOctave);
%             wtb = log10(abs(wtb_temp));
%             Tb = 1./s;
%         else
%                 wtb_temp = cwt(traces(tracein).data,'bump',1/Ts,...
%                     'ExtendSignal',false,...
%                     'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
%                     'VoicesPerOctave',VoicesPerOctave);
%             wtb = wtb + log10(abs(wtb_temp));
%         end
%         count = count + 1;
%     end
    
        p = pcolor(timespec,seconds(T),(wtb));
        caxis([quantile((wtb(:)),0.05),quantile((wtb(:)),1)])
        %p.FaceColor = 'interp';
        p.EdgeColor = 'none';
        colormap jet
        grid on
        hold on
    %     coip = fill([time(1);time;time(end);time(1)],1./[min(s);coib;min(s);min(s)],'w');
    %     coip.FaceAlpha = 0.8;
        ylabel('Period')
        set(gca,'yscale','log')
        grid on
        set(gca,'layer','top');
        %xticks(xticklist)
        xtickangle(0)%, grid minor
        %xticklabels(xticklabellist)
        title(strcat('CWT bump:',...
            datestr(traces(1).startTime,'mm-dd-yyyy')))
        yticks([1:1:5,10:10:50])
    end
    
    %% spectrogram
    % in some cases does bit better than cwt, some worse. could be useful
    % complement to cwt
    %
    if WaveletWeights(3) > 0
        plotcount = plotcount + 1;
        subplot(plotrows,plotcols,plotcount)
    %     count = 0;
    %     specwindow = 120; %seconds
    %     %specnoverlap = 0.9*specwindow; %seconds
    %     specwindow = round(specwindow/Ts/2)*2 + 1;
    %     %specnoverlap = round(specnoverlap/Ts);
    %     specnoverlap = specwindow-1;
    %     for tracein = 1:length(traces) % loop over traces
    %         if count == 0
    %             [sspec,fspec,tspec] = spectrogram(traces(tracein).data,...
    %                 specwindow,specnoverlap,s,1/Ts,'power');
    %             sspec = log10(abs(sspec));
    %         else
    %             sspec = sspec + log10(abs(spectrogram(traces(tracein).data,...
    %                 specwindow,specnoverlap,s,1/Ts,'power')));
    %         end
    %     end
        p = pcolor(timespec,seconds(T),(sspec));
        caxis([quantile((sspec(:)),0.05),quantile((sspec(:)),1)])
        %p.FaceColor = 'interp';
        p.EdgeColor = 'none';
        colormap jet
        grid on
        hold on

        ylabel('Period')
        set(gca,'yscale','log')
        grid on
        set(gca,'layer','top');
        %xticks(xticklist)
        xtickangle(0)%, grid minor
        %xticklabels(xticklabellist)
        title(strcat('Spec'))
        yticks([1:1:5,10:10:50])
        %}
    end
    
    %% reassigned spectrogram
    % doesn't look appriciably different from normal spectrogram
    %{
    subplot(2,3,3)
    count = 0;
    rspecwindow = 120; %seconds
    rspecnoverlap = 0.9*rspecwindow; %seconds
    rspecwindow = round(rspecwindow/Ts/2)*2 + 1;
    rspecnoverlap = round(rspecnoverlap/Ts);
    rs = linspace(1/seconds(Tmax),1/seconds(Tmin),3*traces(1).sampleCount); %evenly spaced freq required for reassigned
    for tracein = 1:length(traces) % loop over traces
        if count == 0
            [srspec,frspec,trspec] = spectrogram(traces(tracein).data,...
                rspecwindow,rspecnoverlap,rs,1/Ts,'power','reassigned');
            srspec = log10(abs(srspec));
        else
            srspec = srspec + log10(abs(spectrogram(traces(tracein).data,...
                rspecwindow,rspecnoverlap,rs,1/Ts,'power','reassigned')));
        end
        count = count + 1;
    end
    p = pcolor(trspec,1./frspec,srspec);
    caxis([quantile(srspec(:),0.05),quantile(srspec(:),1)])
    %p.FaceColor = 'interp';
    p.EdgeColor = 'none';
    colormap jet
    grid on
    hold on

    ylabel('Period')
    xlabel('Time (UTC)')
    set(gca,'yscale','log')
    grid on
    set(gca,'layer','top');
    %xticks(xticklist)
    xtickangle(0)%, grid minor
    %xticklabels(xticklabellist)
    title(strcat('RSpec:',...
        datestr(datetime(traces(1).startTime,'ConvertFrom','datenum'),'mm-dd-yyyy')))
    %}
    
    %% syncrosqueeezed dft
    %just worse resolution and slower that syncrosqueezed wavelet
    %{
    subplot(2,3,4)
    count = 0;
    sstwindow = 120; %seconds
    sstwindow = round(sstwindow/Ts/2)*2 + 1;
    for tracein = 1:length(traces) % loop over traces
        if count == 0
            [ssst,f_sst,tsst] = fsst(padarray(traces(tracein).data,...
                [round(traces(tracein).sampleCount)],'both'),...
                1/Ts, sstwindow);
            ssst = log10(abs(ssst));
        else
            ssst = ssst + log10(abs(fsst(padarray(traces(tracein).data,...
                [round(traces(tracein).sampleCount)],'both'),...
                1/Ts, sstwindow)));
        end
        count = count + 1;
    end
    p = pcolor(tsst,1./f_sst,ssst);
    caxis([quantile(ssst(:),0.05),quantile(ssst(:),1)])
    %p.FaceColor = 'interp';
    p.EdgeColor = 'none';
    colormap jet
    grid on
    hold on

    ylim(seconds([Tmin,Tmax]))
    xlim(round([tsst(end/3),2*tsst(end)/3]))
    ylabel('Period')
    xlabel('Time (UTC)')
    set(gca,'yscale','log')
    grid on
    set(gca,'layer','top');
    %xticks(xticklist)
    xtickangle(0)%, grid minor
    %xticklabels(xticklabellist)
    title(strcat('FSST:',...
        datestr(datetime(traces(1).startTime,'ConvertFrom','datenum'),'mm-dd-yyyy')))
    %}
    
    %% syncrosqueeezed amor wavelet
    %{
    subplot(2,3,4)
    count = 0;
    for tracein = 1:length(traces) % loop over traces
        if count == 0
            [swsst,f_wsst] = wsst(traces(tracein).data,...
                1/Ts, 'amor', 'VoicesPerOctave', VoicesPerOctave);
            swsst = (abs(swsst));
        else
            swsst = swsst + (abs(wsst(traces(tracein).data,...
                1/Ts, 'amor', 'VoicesPerOctave', VoicesPerOctave)));
        end
        count = count + 1;
    end
    swsstplot = swsst;
    swsstplot(swsst == 0) = NaN;
    swsstplot = log10(swsstplot);
    p = pcolor(time,1./f_wsst,swsstplot);
    caxis([quantile(swsstplot(:),0.1),quantile(swsstplot(:),1)])
    %p.FaceColor = 'interp';
    p.EdgeColor = 'none';
    colormap jet
    grid on
    hold on

    ylim(seconds([Tmin,Tmax]))
    ylabel('Period')
    xlabel('Time (UTC)')
    set(gca,'yscale','log')
    grid on
    set(gca,'layer','top');
    %xticks(xticklist)
    xtickangle(0)%, grid minor
    %xticklabels(xticklabellist)
    title(strcat('WSST:',...
        datestr(datetime(traces(1).startTime,'ConvertFrom','datenum'),'mm-dd-yyyy')))
    %}
    
    %% syncrosqueeezed bump wavelet
    %{
    subplot(2,3,5)
    count = 0;
    for tracein = 1:length(traces) % loop over traces
        if count == 0
            [swsstb,f_wsstb] = wsst(traces(tracein).data,...
                1/Ts, 'bump', 'VoicesPerOctave', VoicesPerOctave);
            swsstb = (abs(swsstb));
        else
            swsstb = swsstb + (abs(wsst(traces(tracein).data,...
                1/Ts, 'bump', 'VoicesPerOctave', VoicesPerOctave)));
        end
        count = count + 1;
    end
    swsstbplot = swsstb;
    swsstbplot(swsstb == 0) = NaN;
    swsstbplot = log10(swsstbplot);
    p = pcolor(time,1./f_wsstb,swsstbplot);
    caxis([quantile(swsstbplot(:),0.1),quantile(swsstbplot(:),1)])
    %p.FaceColor = 'interp';
    p.EdgeColor = 'none';
    colormap jet
    grid on
    hold on

    ylim(seconds([Tmin,Tmax]))
    ylabel('Period')
    xlabel('Time (UTC)')
    set(gca,'yscale','log')
    grid on
    set(gca,'layer','top');
    %xticks(xticklist)
    xtickangle(0)%, grid minor
    %xticklabels(xticklabellist)
    title(strcat('WSST:',...
        datestr(datetime(traces(1).startTime,'ConvertFrom','datenum'),'mm-dd-yyyy')))
    %}
    
    %% combined
    %{
    plotcount = plotcount + 1;
    subplot(plotrows,plotcols,plotcount)
%     p = pcolor(timespec,T,(comb));
    A_above_LTA_plot = A_above_LTA;
    A_above_LTA_plot(isinf(A_above_LTA_plot)) = NaN;
    p = pcolor(timespec,seconds(T),A_above_LTA_plot);
    %caxis([quantile((comb(:)),0.1),quantile((comb(:)),1)])
    %p.FaceColor = 'interp';
    p.EdgeColor = 'none';
    colormap jet
    grid on
    hold on
    ylim(seconds([Tmin,Tmax]))
    ylabel('Period (s)')
    set(gca,'yscale','log')
    grid on
    set(gca,'layer','top');
    %xticks(xticklist)
    xtickangle(0)%, grid minor
    %xticklabels(xticklabellist)
    title(strcat('STA-LTA'))
    yticks([1:1:5,10:10:50])
    
    if ~isempty(TQ_tab)
        scatter(TQ_tab_00.start,seconds(TQ_tab_00.T),16,'w','.')
        scatter(TQ_tab_0.start,seconds(TQ_tab_0.T),16,'w','*')
        scatter(TQ_tab_1.start,seconds(TQ_tab_1.T),36,'m','.')
        scatter(TQ_tab_2.start,seconds(TQ_tab_2.T),36,'m','*')
        scatter(TQ_tab.start,seconds(TQ_tab.T),36,'k','o')
%         text(TQ_tab_1.start+10,TQ_tab_1.T,strcat(...
%             'Q=',num2str(TQ_tab_1.Q,'%.1f'),...
%         '  p=',num2str(TQ_tab_1.goodness,'%.2e')),'Color','w')
%         text(TQ_tab_2.start+10,TQ_tab_2.T,strcat(...
%             'Q=',num2str(TQ_tab_2.Q,'%.1f'),...
%         '  p=',num2str(TQ_tab_2.goodness,'%.2e')),'Color','m')
        text(TQ_tab.start+seconds(30),seconds(TQ_tab.T),strcat(...
            'Q=',num2str(TQ_tab.Q,'%.2f')),'Color','k','FontWeight','bold')
    end
    grid on
    %}
    
    %% Wigner-Ville
    % doesn't seem very useful for this purpose
    %{
    subplot(2,3,4)
    count = 0;
    NumFrequencyPoints = 2*ceil(traces(tracein).sampleCount/4);
    NumTimePoints = 2*ceil(traces(tracein).sampleCount/4);
    for tracein = 1:length(traces) % loop over traces
        if count == 0
            [dwv,fwv,twv] = wvd(traces(tracein).data, 1/Ts,...
                'SmoothedPseudo',...
                'NumTimePoints',NumTimePoints,...
                'NumFrequencyPoints',NumFrequencyPoints);
            dwv = log10(abs(dwv));
        else
            dwv = dwv + log10(abs(wvd(traces(tracein).data, 1/Ts,...
                'SmoothedPseudo',...
                'NumTimePoints',NumTimePoints,...
                'NumFrequencyPoints',NumFrequencyPoints)));
        end
        count = count + 1;
    end
    p = pcolor(twv,1./fwv,dwv);
    %p.FaceColor = 'interp';
    p.EdgeColor = 'none';
    colormap jet
    grid on
    hold on
    ylim(seconds([Tmin,Tmax]))

    ylabel('Period')
    xlabel('Time (UTC)')
    set(gca,'yscale','log')
    grid on
    set(gca,'layer','top');
    %xticks(xticklist)
    xtickangle(0)%, grid minor
    %xticklabels(xticklabellist)
    title(strcat('Pseudo W-V:',...
        datestr(datetime(traces(1).startTime,'ConvertFrom','datenum'),'mm-dd-yyyy')))
    %}

    %% hilbert-huang (emd)
    % doesn't seem useful for this purpose
    %{
    subplot(2,3,4)
    count = 0;
    FrequencyResolution = (1/50)-(1/60);
    for tracein = 1:length(traces) % loop over traces
        imf = emd(traces(tracein).data,...
                'Interpolation','spline',...
                'MaxNumIMF',20,...
                'MaxEnergyRatio',100,....
                'SiftRelativeTolerance',0.05);
        if count == 0
            [hshil,fhil,thil] = hht(imf,1/Ts,...
                'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                'FrequencyResolution',FrequencyResolution);
        else
            hshil = hshil + hht(imf,1/Ts,...
                'FrequencyLimits',[1/seconds(Tmax),1/seconds(Tmin)],...
                'FrequencyResolution',FrequencyResolution);
        end
        count = count + 1;
    end

    p = pcolor(thil,1./fhil,hshil);
    %p.FaceColor = 'interp';
    p.EdgeColor = 'none';
    colormap jet
    grid on
    hold on

    ylabel('Period')
    xlabel('Time (UTC)')
    set(gca,'yscale','log')
    grid on
    set(gca,'layer','top');
    %xticks(xticklist)
    xtickangle(0)%, grid minor
    %xticklabels(xticklabellist)
    title(strcat('HHT:',...
        datestr(datetime(traces(1).startTime,'ConvertFrom','datenum'),'mm-dd-yyyy')))
    %}
    
    
    %% waveforms
    %
    plotcount = plotcount + 1;
    subplot(plotrows,plotcols,plotcount)

    %find NPT vert
    stations = {traces.station};
    channels = {traces.channel};
    statchannels = string(stations) + string(channels);
    tracemaxin = find(strcmp(statchannels,'NPTHHZ'),1);
    if isempty(tracemaxin)
        tracemaxin = find(strcmp(statchannels,'NPBHHZ'),1);
    end
    if isempty(tracemaxin)
        %find max trace
        tracemax = zeros(size(traces));
        for tracein = 1:length(traces)
            tracemax(tracein) = var(traces(tracein).data);
        end
        [~, tracemaxin] = max(tracemax);
    end

    plot(time, traces(tracemaxin).data - mean(traces(tracemaxin).data), 'k')
    hold on
%     plot(time, bandpass(traces(tracemaxin).data, [1/120,1/20], 1/seconds(Ts)))
    grid on
    title(strcat('Data:',traces(tracemaxin).network,',',...
        traces(tracemaxin).station,',',...
        traces(tracemaxin).channel))
    ylim([min(traces(tracemaxin).data - mean(traces(tracemaxin).data)),...
        max(traces(tracemaxin).data - mean(traces(tracemaxin).data))])
    ylabel('Velocity (m/s)')
    xlim([min(timespec),max(timespec)])
    %}
    drawnow
    
else
    tqplot = [];

end %end plot

%% compress cwt
if ~isempty(cwt_compress)
    timeedges = timespec(startminind):cwt_compress{1}:timespec(startmaxind+2);
    sedges = 10.^(linspace(log10(1/seconds(Tmax)),log10(1/seconds(Tmin)),cwt_compress{2}));
    spectimebins = discretize(timespec,timeedges);
    sbins = discretize(s,sedges);
    %comb_compress = zeros(length(sedges)-1,length(timeedges)-1);
    [timebinin,sbinin] = meshgrid(spectimebins,sbins);
    timebinin = timebinin(:);
    sbinin = sbinin(:);
    comb = comb(:);
    nanind = isnan(timebinin) | isnan(sbinin);
    timebinin = timebinin(~nanind);
    sbinin = sbinin(~nanind);
    comb = comb(~nanind);
    comb = accumarray([sbinin,timebinin], comb,...
        [length(sedges)-1, length(timeedges)-1], @mean);
    s = (sedges(2:end)+sedges(1:end-1))/2;
    timespec = timeedges(1:end-1) + (timeedges(2:end) - timeedges(1:end-1))/2;
%     if makeplot
%         subplot(2,3,6)
%         p = pcolor(timespec,1./s,comb);
%         caxis([quantile(comb(:),0.1),quantile(comb(:),1)])
%         %p.FaceColor = 'interp';
%         p.EdgeColor = 'none';
%         colormap jet
%         grid on
%         hold on
%         ylim(seconds([Tmin,Tmax]))
%         ylabel('Period')
%         xlabel('Time (UTC)')
%         set(gca,'yscale','log')
%         grid on
%         set(gca,'layer','top');
%         %xticks(xticklist)
%         xtickangle(0)%, grid minor
%         %xticklabels(xticklabellist)
%         title(strcat('compress:',...
%             datestr(datetime(traces(1).startTime,'ConvertFrom','datenum'),'mm-dd-yyyy')))
%     end
end
cwt_struct = struct('wt',comb,'time',timespec,'s',s,'numchannels',count);
    
end %end function

